Matrix product states and variational methods applied to critical quantum field theory 
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We study the second-order quantum phase-transition of massive real scalar field theory with a 
quartic interaction in (1-1-1) dimensions on an infinite spatial lattice using matrix product states 
(MPS). We introduce and apply a naive variational conjugate gradient method, based on the time- 
dependent variational principle (TDVP) for imaginary time, to obtain approximate ground states, 
using a related ansatz for excitations to calculate the particle and soliton masses and to obtain the 
spectral density. We also estimate the central charge using finite-entanglement scaling. Our value 
for the critical parameter agrees well with recent Monte Carlo results, improving on an earlier study 
which used the related DMRG method, verifying that these techniques are well-suited to studying 
critical field systems. We also obtain critical exponents that agree, as expected, with those of the 
transverse Ising model. Additionally, we treat the special case of uniform product states (mean field 
theory) separately, showing that they may be used to investigate non-critical quantum field theories 
under certain conditions. 



INTRODUCTION 

Quantum field theories (QFT) [1] are extremely good 
at describing and predicting the behaviour of fundamen- 
tal particles, as demonstrated by the prediction of the 
Higgs boson over forty years ago and its recent appar- 
ent discovery [2]. Very often, however, obtaining pre- 
dictions from QFT is difficult due, in no small part, to 
the huge Hilbert spaces they are set in. In many cases, 
such as quantum electrodynamics (QED), perturbation 
theory has been used very successfully, yet some impor- 
tant phenomena are not accessible to these methods, no- 
tably confinement in quantum chromodynamics (QCD) 
[3]. Lattice regularizations of QFT's have been very use- 
ful in such cases, often in combination with Monte Carlo 
numerical techniques. Here, however, the sign problem 
[4] presents a challenge - not to mention that such simu- 
lations often require very large computational resources 
to produce useful results (see [5], where hadron masses 
are determined with the help of clusters and a supercom- 
puter) . 

Meanwhile, the numerical study of lattice systems in 
one spatial dimension has benefited greatly from the den- 
sity matrix normalization group (DMRG) method, which 
has been recognized as a variational method producing 
approximate ground states that are matrix product states 
(MPS) [6, 7]. There is a direct link between the di- 
mension of the MPS parameter-space and the amount 
of entanglement that a state can contain, allowing effi- 
cient representation of a great many relevant states [8- 
10], in particular ground states and low- lying excited 
states of gapped systems [11-13], which all lie in the 
low-entanglement "corner" of Hilbert space. Recently, 
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other variational techniques have been applied to MPS 
such as the time-dependent variational principle (TDVP) 
[14], which permits efficient simulation of dynamics, and 
a related excitation ansatz for the determination of dis- 
persion relations [15] for translation-invariant systems in 
the thermodynamic limit. 

Given the great successes of MPS and variational 
methods in studying lattice systems, it is natural to ask 
whether they can be usefully applied to lattice quantum 
fields in (1 -f- 1) dimensions, and whether continuum re- 
sults can be extracted efficiently. A useful test case is 
(/)*-theory which, despite its simplicity, exhibits interest- 
ing behaviour such as spontaneous symmetry-breaking. 
It contains a second-order quantum phase-transition in 
(1 -|- 1) dimensions [16] and is expected to belong to the 
same universality class as the transverse Ising model [17] 
so that critical exponents should be the same for both. 
In fact, DMRG has already shown promise when applied 
to (^'^-theory [18, 19], reproducing the expected critical 
behaviour and obtaining values of the critical parameter 
close to those of Monte Garlo studies [20-22] . 

We use variational methods with MPS to obtain the 
ground-state field expectation value and low-lying exci- 
tation energies of 0^-theory. The scaling of these quan- 
tities in parameter-space allows us to locate the critical 
point and to determine the critical exponents. We be- 
gin by introducing QFT and real scalar field theory in 
section I, showing how it can be put on a spatial lattice 
and discussing its critical behaviour. In section II, we 
define the uMPS variational class and the correspond- 
ing TDVP algorithm and excitation ansatz before detail- 
ing our variational conjugate-gradient method for finding 
ground states. Section HI is the main part of this work, 
in which we apply these techniques to (/)^-theory and ob- 
tain our estimate for the continuum critical parameter, 
which we compare with previous results from the liter- 
ature, as well as values for critical exponents and the 
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central charge. We also separately assess the usefulness 
of mean-field theory (a special case of MPS) for study- 
ing QFT, which is an attractive tool because of the low 
computational complexity required to estimate physical 
quantities. 

We have kept the software developed for this work in- 
tentionally general such that it may be of use to others. 
It is available under a permissive open-source license [23] . 



I. QUANTUM FIELD THEORY 

We introduce the basic principles of quantum field the- 
ory using the Hamiltonian formulation with real scalar 
fields as an example, defining interacting (^"^-theory both 
in the continuum and on a spatial lattice. We then dis- 
cuss its spontaneous symmetry-breaking, which corre- 
sponds to a second-order quantum phase-transition in 
(1-1-1) dimensions. The beginning of this section is 
based partly on lectures on quantum field theory given 
by Marco Zagermann at Leibniz Universitat Hannover in 
2010/11 and also on [1]. 



A. Real scalar field without interactions 

Quantum fields in Minkowski space-time are quan- 
tum systems set in an uncountably infinite-dimensional 
Hilbert space which can usually be divided naturally into 
subsystems corresponding to points in momentum space. 
They can often be constructed from a corresponding clas- 
sical field defined by a Lorentz-invariant action. The clas- 
sical field is then quantized in such a way as to produce a 
consistent Hilbert space and Hamiltonian where Lorentz- 
invariance and causality are maintained. 

As an example, take a classical real scalar field 4>{x) G 
M with action 



(1) 



where v = Q . . . d — 1, d,j = djdx^ and scalar products 
are defined via the Minkowski metric with the "mostly- 
minus" signature (1,— I-- - — 1). We use x to denote 
the spatial part of a Minkowski vector x. The integrand 
£[0(a;), 9(/'(a;)] is called the Lagrangian density. Stability 
with respect to a small variation (5(/) requires that the 
Euler-Lagrange equations 



dC 







are satisfied, leading in this case to the equation of mo- 
tion 

(□ + /io)'^(a^) =0, 

where □ — d^d'^, which is the Klein-Gordon equation. 
The lack of non-linear field-terms in the equations of mo- 
tion makes this a free (non-interacting) field. Performing 



a Fourier transform (j){t,x) = J ^2^^1-1 e'P-"'0(t, p), we 
can rewrite the equations of motion as 

which has the form of a simple harmo nic oscil lator with 
angular frequency w(p) = E{p) = \/ + Mo- '^^^ 
thus think of the classical Klein- Gordon field as a set of 
independent harmonic oscillators, one for each point in 
momentum space. 

To quantize this free scalar field theory, we start from 
the classical Hamiltonian. The canonical conjugate mo- 
mentum 7r(x) corresponding to the coordinate <i>{x) is 



'k{x') 



dL 



d{do^{x)) 



= do<i){x) = (j){x) 



and, performing a Legendre transformation, the Hamil- 
tonian density is 

The coordinate 4>{x) and the momentum 7r(x) obey the 
Poisson-bracket relationship 

{(t){t,x),TT{t,y)} ^ S{x ~ y), 

where S{x) is the {d — l)-dimensional Dirac delta distri- 
bution. The ingredients required for a canonical quanti- 
zation of the classical theory are now ready. To proceed, 
we replace the classical phase-space coordinates in the 
above relations with operators (one for each space-time 
coordinate) obeying the commutation relation 

[<j){t,x),n{t,y)] ^ i6{x - y). 

The field operator (f) = (f)^ \s Hermitian because the clas- 
sical field was real-valued. We now have an operator- 
valued field where the operators 4>{x) must obey the 
Klein-Gordon equation. The harmonic oscillator picture 
of the classical field suggests attempting to write solu- 
tions in terms of quantum harmonic oscillators. Making 
the same move to momentum space as before, we can 
write a general solution as a superposition of plane waves 
using Fock-space creation and annihilation operators 
and Or, 



dp 



\pO=E(p), 



where [ap,a^] = (27r)'* ^5{p — q) and p° = E{p) — 

\/ p^ + Mo- The Hilbert space contains the vacuum 
flp |0) — Vp and countably infinite excited states 
(a^)" |0) for each momentum p S E"*"^. Using 7r(x) = 

(bix). we can write the Hamiltonian as 



U = 



dp 



(27r)« 



1, 
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where the second term does not annihilate the vacuum, 
leading to an infinite vacuum energy contribution. This 
is perhaps not too surprising: We are summing up an infi- 
nite number of ground state energies, one for each Fourier 
mode, each of which is the energy contained within an 
infinite volume of space. Since it is energy differences 
that are observable, and not absolute energies, this in- 
finite contribution should not cause any problems. A 
general eigenstate a|,a^ ■ ■ • |0) has energy (ignoring the 
infinite vacuum contribution) E{p) + E{q) + . . . and is 
also a momentum eigenstate with momentum p + </ + ... , 
where the momentum operator can be obtained via the 
classical theory as the conserved quantity associated with 
spatial translations (using Noether's theorem) . The field 
operator (t>{t, x) (in the Heisenberg picture) acts on the 
vacuum to create a superposition of momentum eigen- 
states resulting in a particle localized at the space-time 
coordinate x. 

A quantity that turns out to be very useful is the two- 
point correlation function {Q\(j){x)(j){y)\Q) , which can be 
interpreted as the probability of a particle created at 
point X propagating to point y (or vice-versa, depend- 
ing on the time-coordinates). For this reason it is also 
called the propagator. It has the form 



D{x ^ V) ^ {0\cp{x)c^{vm 
dp 1 



(27r)rf-i 2E{p) 



A related quantity is the Feynman propagator, which is 
defined as 



Df{x - y) 



dp 



(O|T0(a;)0(y)|O) 



' D{x — y) for x'^ > 



(2) 



D(y ~ x) for a;° < y° ' 



where T denotes the time-ordered product and the rela- 
tion to Dix — y) can be found using contour integration, 
with the infinitesimal shift e providing a prescription for 
treating the poles. Dp is a Green's function of the Klein- 
Gordon equation 

(□ + fil)DF{x - y) = -iSix - y). 

The integrand oi Dp has poles given by the mass pa- 
rameter at p2 = Mo- For an interacting theory, the poles 
no longer correspond to the mass parameter /Iq, but are 
shifted away from this point due to self-interaction. The 
shifted poles of the propagator then correspond to the 
physical mass of a particle whereas the "bare" parameter 
Mg does not. That the poles of the propagator correspond 
to the particle mass can be seen by inserting the iden- 
tity, written in terms of the (unspecified) eigenstates of 
the Hamiltonian (interacting or not) and the momentum 
operator, into the expression for the propagator. The 
identity thus formed is 



dp 



{2tiY-^ 2E{p,X) 



|Ap) (A, 



where \fl) is the vacuum, |Ap) is the zero-momentum 
state |Ao) boosted to momentum p, and E{p,X) = 



with TO^ being the mass or energy of |Ao) 



Evaluating {n\(l>{x)I(l){y)\fl) leads to the Kallen-Lehmann 
spectral representation of the Feynman propagator (see 
section 7 of [1]) 



Df(x - y) 



°° dAP 



27T 



Po{M')DF{x-y,M'), 



where Dp{x — y,M'^) is the Feynman propagator with 
mass-parameter A'P (instead of /x§) and 

Pp(M^) ^ ^(2^)5(M2 - ml)\ mmK) ? (3) 



is the spectral density. We drop the (fi| term, since 
it adds at most a constant term to the propagator. 
Given that the theory contains single-particle states. 



the spectral density contains a pole at M 



where 



m is the mass of a single particle, followed by a gap be- 
fore further excitations appear. In this case, the Feynman 
propagator can be separated into a one-particle contribu- 
tion and the rest. With a Fourier transform we have 



da; e'P-^Di.(.T-y) 



iZ 



— mn? 



where Z is a real number coming from the | (ri|(/)(0)|Ao) P 
factors. The single-particle term has a pole at — m? , 
with the other terms showing up at higher momenta. 



B. Interacting fields 

So far, we have considered the quantized free scalar 
field, whose solutions are plane-waves. A free (non- 
interacting) field is, however, not directly relevant to 
physics, since a lack of coupling implies a lack of mea- 
surable consequences. We can introduce interactions by 
adding a term to the Lagrangian density C — C^cc + >Cint 
that leads to non-linear equations of motion. An exam- 
ple for real scalar field theory, and the case we focus on 
in this paper, is the quartic interaction term 

where A is the coupling constant, or the strength of the 
interaction. The resulting theory is often referred to sim- 
ply as "(/)'*-theory" . Its equation of motion is 

{u + ^,l)d,{x)^-J^<j,^ 

which no longer has simple plane-wave solutions. Since 
interacting theories are generally difficult or impossible 
to solve analytically, other approaches such as discrete 
(lattice-based) numerical simulation or perturbation the- 
ory are needed. The perturbative approach is used to 
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study scattering, where it is assumed that the incom- 
ing and outgoing states far from the scattering location, 
called asymptotic states, can be described by the non- 
interacting field theory. Scattering is then represented 
by a unitary operator, the "S-matrix" , relating the in- 
coming and outgoing states. Elements of the S-matrix 
can be calculated perturbatively in powers of the cou- 
pling constant. The individual terms in the expansion 
have a regular form and can be conveniently represented 
using Feynman diagrams. 



Since, for this work, we perform numerical simulations 
on a lattice, we do not go into perturbative calculations in 
detail. As mentioned above, the perturbative calculation 
of the propagator in an interacting theory reveals a shift 
of the pole mass away from the bare mass parameter /Xq, 
resulting in a different, "dressed" physical mass Mphys- 
The mass shift is due to the interaction of the field with 
itself, which can involve modes of any momentum. In 
fact, taking all possible momenta into account, the shift 
diverges. Introducing a momentum cut-off into calcula- 
tions, for example via a lattice, makes the shift dependent 
on this cut-off. Since the physical mass of particles can- 
not diverge, and because the bare parameter is not itself 
measurable, the bare mass is adjusted such that the pole 
of the propagator has the correct (measured) value, even 
if this means that the bare mass diverges. The procedure 
of adjusting bare parameters to cancel contributions from 
self-interaction is called "renormalization" . In general, 
the bare mass is not the only parameter that must be 
renormalized. Others, such as coupling constants, may 
also be affected. 



The need for renormalization and the presence of diver- 
gent shifts can be interpreted as signs that the theory in 
question is an effective low-energy limit of a more funda- 
mental one [1]. The momentum scale where the effective 
theory breaks down then becomes a natural cut-off, such 
that divergent quantities are avoided. In the standard 
model of particle physics, a candidate for this cut-off is 
the Planck scale ^ 10^^ GeV, where gravitational effects 
are expected to play a significant role. However, since we 
don't know which theory describes physics beyond the 
standard model, we also cannot know the exact location 
of the cut-off, which may occur at far lower energies. 



to first order in A is 



.A 



- Mo 



- Mo 



dq 



{2nYq^-^li + ie 



0(A2), 



(4) 




where the first term is the free-field propagator (2) and 
the second term is the first-order correction. The part in 
square brackets —\5^\ diverges. 

The name "one-loop" 
for the correction term q 
comes from the cor- 
responding Feynman 
diagram (Figure 1), where 

the two free-particle prop- ^ V ^ ^ 

agator factors outside > — ^~ 

the square brackets cor- 

respond to incoming and ^^6^"^^ Feynman diagram 
, . , . 1 , of the one-loop correction to 

outgomg particles, each ^, ^ ^. , 

. or , . , the tree particle propagator 

with momentum p, which ^^-theory 
are represented by incom- 
ing and outgoing lines in the diagram. The integral over 
q inside the brackets corresponds to a "virtual" particle 
and is represented by a loop. Additional loop-terms 
appear at higher orders. Another way of writing the 
propagator, recognizing that the perturbative expansion 
in further loop terms results in a geometric series, is as 



dxe''P''DF{x-y) 



Mo 



where the mass shift S^^ now contains all the loop- 
corrections. To first order, S^^ = and the shift di- 
verges. Higher-order contributions to do not diverge 
in (1 -|- 1) dimensions, such that removing the divergence 
is already achieved by adjusting the bare parameter /ig by 
A finite shift coming from higher-order corrections 
remains, but is unimportant for the purposes of investi- 
gating critical behaviour, where the physical parameters 
must merely be well-defined. There are also finite renor- 
malization factors corresponding to the field and the 
coupling A, which we ignore for the same reasons. 



C. Real scalar field theory on a lattice 



For our purposes, it is sufficient to briefly examine the 
only divergent (in the absence of a cut-off) term in 
theory in (1 -I- 1) dimensions [24], which is the "one-loop" 
correction to the propagator. The propagator describes a 
simple "scattering" event involving a single incoming and 
outgoing particle, which we can examine using the same 
perturbative methods as are used for more complicated 
scattering events. The Fourier-transformed propagator 



Discretizing the space in which a quantum field lives is 
a possible way of making field theories accessible to non- 
perturbative methods. It corresponds to a dramatic re- 
duction in the dimension of the Hilbert space and implies 
a momentum cut-off, making loop-integral contributions, 
which may be divergent in the continuum, finite on the 
lattice. 

We use a spatial discretization to permit the use of 
matrix product states (MPS). Time remains continuous. 
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as we simulate dynamics using the Hamiltonian formal- 
ism. For more details of this procedure, see [25]. We 
also work with an infinite lattice (in the thermodynamic 
limit) since this is possible using MPS and is a more 
realistic setting for a field than a finite lattice. By as- 
suming spatial uniformity of ground states, the number 
of variational parameters needed to approximate states 
remains manageable. Our lattice version of the classical 
continuum theory introduced in (1) is given, in (1 + 1) 
dimensions, by the Lagrangian 



L = a 



E 



Mo -2 



2a2 



where the sum is over the lattice sites, a is the lattice 
spacing, and the spatial part of the derivative term has 
been replaced by a finite difference. Letting a —>■ re- 
covers the Lagrangian of the continuum free scalar field 
theory. Applying the Euler-Lagrange equations results 
in 

I- ^(20„ - <j)n-l " 4>n+l) + fJ-l4>n = 0, 



where the term in brackets can be interpreted as the sec- 
ond derivative on the spatial lattice. As with the classical 
continuum theory, a Fourier transform 



dp° 



IT / a 



TT / a 



27r 



ip X ip na 



diagonalizes the equation of motion: 
dpO r'" dpi 



2lT 



— tt/ a 



2ti 



iP 



0\2 



Mo 



0(p) = 0. 



The Fourier-transformed Green's function G{p, a) satis- 
fies 



so that 



G(p,a) = 



2 

Mo 



G{p,a) 



(p0)2 

which, in the limit a 



■A sin^ 



a 



Mo 



0, becomes 
i 



Gip) = 



9 ' 

Mo 



in agreement with (2). 

Moving now to the quantized and interacting t/)^- 
theory, we can write down the one-loop correction to the 
physical mass by analogy with (4) 

- iS^l = 

.A r dp" r'" dpi i 
"'2 ./ '2^ . L^u 



— -K I a 



(p0)2-^sin2(^^) 



Mo 



which again agrees with (4) as a — > 0. Integrating over 
p^ using contour integration, this becomes 

A dpi 1 



A r''^ 

4<5^? = -i-y 

^ J —7T fa 



27r 



M§ 



which can be written in terms of the complete elliptic 
integral of the first kind 



r/2 



1 



K{k) = / d6- 

'0 «/l - /fc2sin2(6l) 



This leaves 



-iSiil 



.A 1 




(5) 



which is convenient for calculation using numerical com- 
puting packages, where the elliptic integrals are com- 
monly implemented as high-accuracy approximations. 

To investigate behaviour using the time-dependent 
variational principle, we need the Hamiltonian form of 
the interacting lattice theory. With the interaction term 



and using 7r„ 



is 



2a? 



dL 

90„ 



2a2 



the Hamiltonian 



Mo 



A 

4!' 



The parameters A and /^q have dimension [mass] 2. Re- 
placing them with dimensionless quantities /ig = ij^qO^ 
and A = Aa2 allows us to write the dimensionless Hamil- 
tonian H = Ha as 



+ 



Mo 



A 

4!' 



eliminating the explicit appearance of a. Adjusting the 
lattice^ spacing now corresponds to altering the param- 
eters A and ji^. This dimensionless form is convenient 
for finding the continuum limit of the quantum critical 
theory (see section HI). 

Noting that H takes the form of a many-body Hamil- 
tonian with a nearest-neighbour interaction, it becomes 
natural, especially with regard to the later use of ma- 
trix product states, to use a basis given by position- 
space creation and annihilation operators [a„, a\^] = 5nrm 
On |0) = to define the quantized lattice theory. The field 
operator and the conjugate momentum operator can then 
be defined as 



V2' 



[al 



an) 



and 



V2 



such that the desired equal-time (Schrodinger picture) 
commutation relation 



, 71'm] = i<5n 



is satisfied. 
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D. Spontaneous symmetry- breaking 



The (/(^-theory action 
1 



dx 



4! 



(6) 



is manifestly invariant under the discrete transformation 
(/) —(f>. A given state may or may not share this sym- 
metry. Should the ground state of a QFT break a sym- 
metry of the action for some set of parameters, the theory 
is said to exhibit spontaneous symmetry-breaking. The 
word "spontaneous" refers to the fact that there are then 
multiple ground states (the number of ground states is 
equal to the order of the symmetry), such that the ac- 
tual ground state of the system, obtained for example by 
cooling, makes a seemingly spontaneous "choice" . 

Classically, the ground state lies at the minimum of a 
potential function. The (f>*-ih.eoiy action (6) contains the 
classical effective potential 



A 

4!' 



which, for /Xq > 0, has a single minimum at = 0, 
leaving the symmetry </> — ^ —cj) intact. However, with 
/i§ < there are two minima and hence two distinct 
ground states at ±0o,ci > that break the symmetry, as 
illustrated in Figure 2. 




Figure 2. The classical effective potential in ^'^-theory illus- 
trated for /io > (red) and fj,Q < (blue) showing the two 
possible ground states for the latter case. 

The symmetry-breaking persists in the quantized the- 
ory in (1 -|- 1) dimensions, which possesses symmetric 
and symmetry-broken phases distinguished by the vac- 
uum (ground state) expectation value of the field oper- 
ator (f2|(/)|ri), henceforth abbreviated to {(p). Since the 
bare mass parameter /Iq must diverge in the continuum 
in order to renormalize the physical mass (see section 
IB), the relevant parameter is not /Zg as in the classical 
case, but the renormalized mass /z^, where we use the 
definition 

2 2,12 

with S^i being the one- loop correction defined in (4). 
Note that //^ is distinct from the physical mass {fij^ ^ 
A'phys) additional finite corrections. Since we are 



not generally working under weak coupling conditions, 
the usual perturbative calculation of the full mass shift, 
and hence the physical mass, is not applicable, so that the 
physical mass cannot be used as a parameter (although 
it can be found numerically on the lattice). 

When moving through parameter space (A, /i|j), the 
transition between the asymmetric ground-state phase 
and the symmetric phase represents a second order quan- 
tum phase transition [16] with order parameter (0). 
There therefore exist critical points in (A, /i^) at which 
the theory becomes massless and scale-invariant (the 
correlation- length ^ becomes infinite). A scale- invariant 
theory should be described by dimensionless parameters, 
yet the two parameters A and fij^ have dimension [mass]^ 
in (1 -|- 1) dimensions. As such, the proper parameter 
must be the ratio A//x|j. This means there is a line in 
parameter-space corresponding to the critical theory. 

The lattice theory also contains critical points, where 
the critical parameter A//i|j^ = ^IP^R c ^'^^ depends on 
the lattice spacing a (which defines a momentum cut-off) 
or, equivalently, on A, so that we may write A//i|.g(A). 
This dependency is expected to be logarithmic due to 
infrared corrections in the critical theory [26] where the 
physical mass goes to zero. Such a dependency has been 
observed in Monte-Carlo simulations [21]. To obtain the 
critical parameters of the continuum theory, we take the 
limit of X/jl\ j,(A) as A — >■ 0. 

Note that it is the lattice correlation length ^ — £,a^^ 
that goes to infinity at the critical points of the lattice 
theory. For this reason, there are two possible interpre- 
tations of the lattice critical point: Either as a lattice ap- 
proximation a > to the continuum critical point where 
^ — cx), or as a continuum limit a — > of a non-critical 
theory ^ < cx). Since we are interested in the critical 
continuum theory, we will always use the former inter- 
pretation. 

In the vicinity of the critical point, physical quanti- 
ties scale according to power laws (see [27] or section 13 
of [1]). For the order-parameter (</>), in the symmetry- 
broken phase where ((/>) 7^ 0, we can thus expect 



(<^) = A{\) 



A 



A 



/3(A) 



where A{\) is some constant and /3(A) is the critical ex- 
ponent. We also define the scaling for the energy (or 
mass) of the lowest-lying excitation; 



Ai; = B{\) 



A 



The energy AE should correspond to the particle mass 

/^phys 

(given by poles in the propagator) in the symmet- 
ric phase, but may belong to a topologically non-trivial 
soliton (kink) excitation in the symmetry-broken phase 
(providing a localized transition — — between two 
different ground states at x ±00). 
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Predictions for the critical exponents /3 and v in the 
hmit A — >■ can be obtained based on the universahty 
principle, which conies from renormalization group the- 
ory [27]. (/)'*-theory in (1 + 1) dimensions has been shown 
to be a continuum limit of the transverse Ising model [17] 
and, as such, is predicted to share its critical exponents. 
For (0) we thus expect /3 = 1/8 and for Ai? we expect 
V = \. For more information on the critical behaviour 
of (^^-theory, including a derivation of these critical ex- 
ponents, see [28]. The critical parameter A//i^ j, is not 
a universal quantity, depending instead on the particu- 
lars of (^''-theory. It is also not accessible to perturbative 
techniques [16], making it an interesting target for lattice 
methods. We estimate it, as well as the critical exponents 
defined above, in section III. 



II. MATRIX PRODUCT STATES 

Matrix product states (MPS) are pure states of one- 
dimensional lattice systems with a particular form that 
puts a limit on the amount of entanglement a state can 
contain. The amount of entanglement is related to the 
bond-dimension Z?, which is the dimension of the matri- 
ces that make up the state coefficients. Most quantities 
(assuming open boundary conditions) can be calculated 
with complexity 0{ND^), where N is the number of lat- 
tice sites or, in the case of uniform (translation invariant) 
MPS in the thermodynamic limit (uMPS), the number 
of necessary solver iterations. 

In this section, we define uMPS and derive an im- 
plementation of the time-dependent variational princi- 
ple (TDVP) as well as a related ansatz for determin- 
ing excitation energies. Both algorithms were first de- 
scribed by Haegeman et al. [14, 15]. We also set out a 
variational conjugate-gradient algorithm, demonstrating 
significantly improved convergence speeds for ^'^-theory 
compared to the TDVP. 



A. Uniform MPS in the thermodynamic limit 



Uniform matrix product states (uMPS) have the form 



tions 



\^{A)) - 5] 4 

{s}=0 



+ 00 



n 



where |s) = |. . . si . . . SAf . . .) and the site-independent 
d X D X D tensor A contains the parameters for the en- 
tire state. The boundary vectors vl and vr are of length 
D and are irrelevant in calculations so that we may ignore 
them. The uMPS states form a sub-manifold of Hilbert 
space A^uMPS C H that depends on D. They have dD"^ 
complex parameters, but possess one non-physical degree 
of freedom corresponding to the norm and — 1 gauge 
degrees of freedom due to invariance under transforma- 



(7) 



where the trivial transformation g = clis not counted, so 
that the number of physical degrees of freedom is {dD^ — 
1) — (Z?^ — 1) = D^id ~ 1). The norm is determined by 
the infinite power of the x matrix E = ^^A'^i^A'' 
so that the spectral radius of E must be one: p{E) — 1. 
We further require that E has a unique eigenvalue of 
greatest magnitude, which must then be equal to one in 
order to obtain well-defined expectation values and to 
avoid dependencies on the boundary vectors [29]. 

The left and right eigenvectors of E with eigenvalue 
1 we name {l\ and \r) respectively. Via the Choi- 
Jamiolkowsky isomorphism, we may also define D x 
D matrices / and r such that A'^HA^ ~ I and 
A^rA*^ = r. Numerical computation of quantities 
involving E is more efficient in this matrix representa- 
tion, scaling with 0{D^) rather than 0{D^) (assuming 
a naive matrix- multiplication algorithm). Single-site ex- 
pectation values can be computed as 



{^{A)\o\^{A)) = {l\E''\r)^tT 



lJ2^'rA'^ {s\o\t) 



where E° = J2s,t (sW) A* (g) A''. 

The gauge freedom (7) implies that there is no unique 
MPS representation of a given state. There are, however, 
useful forms for the uMPS tensor A, of which the so- 
called right canonical form has the properties 



A' A''' =Id r = lD and 

s 

lal3 = SapXa (a, (3 = 1 . ■ ■ D) , 



(8) 



where Aq, are the Schmidt coefficients corresponding to 
decomposing the system into two infinite halves j^*) = 
J2a=i '^a 'S^ \ipR), with orthonormal Schmidt vectors 
for the left and right halves and \iPr). This makes 
explicit the relationship between the amount of entangle- 
ment possessed by a uMPS state and the bond dimension 
D, which is equal to the Schmidt rank of the half-chain 
decomposition. Since the Schmidt coefficients are also 
the eigenvalues of the density matrix corresponding to 
the reduced state on the half-chain, we can easily calcu- 
late the corresponding von Neumann entropy as 



S 



(9) 



a=l 



The conditions (8) fix all gauge degrees of freedom and A 
can always be made to fulfil them by performing a gauge- 
transformation. This can be verified using the eigenvalue 
equations E\r) = \r) and {l\ E — We find that a 
gauge transformation affects I and r as / — > g^^Ug^^ 
and r grg^ which, together with (8), fully specify the 
g needed to put an arbitrary A into canonical form. 
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To implement the TDVP for uMPS, we need to un- 
derstand the tangent plane T|,j,(yi)^ to A^uMPS at a point 
1 4" (A)). Uniform tangent vectors have the form 



\HB))=Y,B,\d,^>{A)) 



(10) 



+ CX3 d—1 
n=-oo {s}=0 



n 



+00 



n 



where we use the shorthand notation \di^{A)) = 
d/dAi \^{A)) with the index i running over all entries 
of the tensors A and B. We call the tensor B the 
parameter-space tangent vector. Changing the state pa- 
rameters as A — )■ A-\-dTB changes the state as |^'(A)) ^■ 
\^{A)) + dT|$(B)). We also define "boosted" tangent 
vectors for uniform systems 



+ 00 



d-1 



E E -I 

n=-oo {s}=0 



n 



B' 



+ 00 



n 



_i—n-^l 



(11) 

vr \s) , 



representing different momentum sectors p such that 
|$o(5)) = These are useful for studying exci- 

tations (see section II C). 

As with the state |^'(^)), there are non-physical de- 
grees of freedom in the parameter tensor B. Apart 
from the state itself lying in the tangent place |^(A)) G 
they also include infinitesimal gauge transfor- 
mations such that a tangent vector is invariant under 
|$p(B)) ^ \%{B+Up{x))} with J\f^{x) = e-^PxA' - 
A'^x. That \^p{J\fp{x))) corresponds to an infinitesi- 
mal gauge transformation can be checked by using one- 
parameter site-dependent gauge transformation matri- 
ces gniv) = I + jyxe'P" to transform the state {Af^ — > 
gn-iAj^g~^), taking the derivative d/dr]\^{A))\^^Q to 
obtain the infinitesimally transformed state. 

All non-physical degrees of freedom can be eliminated 
by requiring that B satisfy a gauge-fixing condition such 
as the right gauge-fixing condition 



(12) 



Note that ($p(B)|«'(A)) = 2tt6{p) {1\E^\t) where the 
second factor is zero according to (12). Hence, this con- 
dition also includes orthogonality to the ground state 
|\E'(A)) for momentum zero, which cannot be obtained 
by a mere gauge transformation. If we start with an ar- 
bitrary B, then for momentum zero we have to manually 



impose {l\Ej 



(orthogonality to the ground state) , 



after which we can bring it into a form where it satisfies 
(12) by doing a gauge transformation. For non-zero mo- 
mentum, a gauge transformation alone is sufficient. We 
can see this by making the replacement B ^ B +J\fp{x) 
in (12), resulting in 

\xr) = {E-le-'Pr'Ef\r). 




The time-dependent variational principle 

We wish to compute ., 
the time evolution of a 
quantum state. Because 
the dimension of the 
Hilbert space is large, 
we restrict ourselves to 
a class of relevant states 
I 'I' [a]) with a manage- 
able number of param- 
eters a e C'', d <C dim('H). This defines a sub- 
manifold Men. Given a starting state |5'[a(t)]), the 
Schrodinger equation gives us the infinitesimal evolution 
\'^{t + dt)) = l'I'[a(t)]) ~ iH\<l'[a{t)]), where generally 
\<l'{t + dt)) ^ M because the step -\H \<i[a{t)\) (blue 
arrow) need not lie within the tangent plane T to M. To 
optimally approximate time evolution whilst remaining in 
M, we project the exact step onto T, which means find- 
ing a tangent vector |$) £ T (red arrow) that minimizes 
\\iH \^[a{t)]) + 1$) Ip. A tangent vector has the form 
Mb]) = b> (with = d/da^ !*[«])), leading to 

the flow equations 



i|$[a(t)]) = i9,*)ff^''-{9,*lH!*) 



where g-' is the inverse of the puUback metric Qj^ ~ 
{dj'i'ldk'i') (assuming gjk has no kernel). We identify 
Idj'i/) g-'^ {dk'i] as the projector onto T. As a simplifi- 
cation, we have taken \'^{a)) to be always normalized 
and to be a holomorphic function of a. 



Figure 3. 

which we can solve to obtain x. For p ^ 0, the solution is 
unique (assuming r is full-rank). In casep — the inverse 
must become a pseudo-inverse, leaving freedom x ^ x + 
cl corresponding to the null space of {E — T), which we 
eliminated from E^ \r) by imposing orthogonality to the 
ground state. However, for p = this freedom in x is not 
part of the gauge group: Afo{cl) = so that the condition 
fixes exactly the gauge (and norm for p = 0) degrees of 
freedom. Restricting B so that it always satisfies (12) 
can be achieved using the parametrization 



B'{x) = r^/^xF^-i/^^ 



(13) 



where x € Af£)x_D(!i-i) and the D{d — 1) x dD matrix 
[V]ia.s):p ~ [V'^]ap (whcrc the index (a,s) combines s 
and a) is defined so that contains an orthonormal 
basis (VV^ = I) for the null-space of R\ with 



resulting in VR — 0. 



B. Time-dependent variational principle for matrix 
product states 

To apply the time-dependent variational principle 
(TDVP — see Figure 3) to uMPS we have to find an 



9 



X that satisfies 

X = argmin|||$(B(a;'))) +iH\'i'{A))\\ 

x' 

for a given where B(x) is the gauge-fixing 

parametrization (13) and |$(_B)) is a uniform tangent 
vector as defined in (10). We minimize the expression by 
setting its derivative with respect to x^ equal to zero. To 
do this, we need to calculate the two terms containing 
x^ . The first is the tangent vector norm 

i^={^{B)\'i>{B))=WB'g,j, 

where the summation indices i and j run over all the 
entries of B. With the gauge-fixing parametrization, this 
simplifies to 

T]{x) = {<^{B{x))MB{x))) = |Z| tr [x^x] , (14) 

where |Z| represents the size of the infinite lattice. We 
also have the Hamiltonian term 

{HB)\H - {H) |*(A)) = W (d,,m - (H) \^{A)) , 

where we are free to subtract (H) = {^{A)\H\'ii{A)) 
without changing the result due to (12), which ensures 
(<f>(i?)|\l/(A)) = 0. Assuming the Hamiltonian is uniform 
and can be written as a sum of nearest-neighbour terms 
H = J2n hn,n+i, this simplifies to 

mB{x))\H-{H)\^{A)) = \Z\tv [x^F], 

with 

S 

K contains the sum of Hamiltonian terms over one half 
of the infinite lattice 

-)-oo 

\K) = Y,{ErE'i^\r), 

with 

u,v 

so that 

^AB = C''*(8)An3t <^ 



represents a single term in — (H) acting on a pair of 
sites. Since E has a unique eigenvalue of largest magni- 
tude with value 1, we can split such infinite sums into 
two parts 

+ C30 

n=0 

with the projector Q = Q" = I — \r) {l\ leading to 
p{QEQ) < 1, turning the first term into a geometric 
series 

J2 QiQEQTQ - Q(I - QEQ)-^Q, 

n=0 

whilst the second term is zero due to {l\E^^\r) = 
(*(A)|/i - (h) |*(A)) = 0. We thus have 

\K) = Q{1 - QEQ)-'QE'^A \r) = (I - EfE^^ \r) , 

where P denotes the pseudo-inverse. \K) can be calcu- 
lated directly, but would involve 0{D^) operations. In- 
stead, we avoid the inverse by re-arranging to give 

(I - QEQ) \K) = QE'Xa \r) , (15) 

which can be solved in the matrix representation for K 
with complexity 0{D^) using a sparse solver. 
Finally, we obtain the TDVP flow equations 

A' = -iB'{F), 

giving us the time-evolution of |^'(^)) € A^uMPS that 
best approximates the exact (Schrodinger) evolution. We 
may integrate them numerically using the Euler method 
with the following algorithm: 

1. Calculate F (including prerequisites C, K). 

2. Take a step by setting A^it + dt) = A^{t) - 
idtB^iF). 

3. Restore canonical form of A using a gauge trans- 
formation. 

4. Compute I and r and normalize, then compute 
other desired quantities, such as the energy, and 
adjust the step size dt as required. 

Normalization is necessary despite gauge-fixing because 
we take finite time steps along tangent vectors. For the 
same reason, the gauge degrees of freedom will also drift 
so that we must perform a gauge transformation if we 
wish to maintain canonical form (which reduces compu- 
tational requirements due to the simple forms of / and 
r). Determining the eigenvectors {l\ and \r) can be done 
iteratively (using a sparse eigensolver) with per-iteration 
complexity 0{D^). If the corresponding eigenvalue is 
not 1 then A should be scaled appropriately to normal- 
ize the state. The total complexity of the algorithm is 
0{nit,:D^), where nitr is the number of iterations required 
to find (/| and |r) plus the solver iterations needed to ob- 
tain K using (15). 
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1. Imaginary time evolution 

Imaginary-time evolution can be seen as a gradient- 
following minimization method applied to the energy 
functional H{^,'^) = Taking the first deriva- 

tive with respect to (^'| results in ^ H\^) so that 

a small step dr in the direction —H\'^) should take us 
closer to the ground state (given that \^) is not orthogo- 
nal to it). The same result is obtained by replacing t with 
—IT in the Schrodinger equation, hence "imaginary-time 
evolution" . 

The TDVP flow equations can be used to efficiently ap- 
proximate the exact imaginary-time evolution by making 
the same replacement t — —it. If we start with a state 
in some variational class that is not orthogonal to the 
exact ground state, integrating the flow equations will 
then locate the best ground-state approximation within 
the class unless we get stuck in a local minimum. The 
norm 77 (defined in (14)) of the approximate evolution 
vector 1$) acts as a convergence measure: It represents 
the size of the gradient H\^) as projected onto T, which 
goes to zero at the energetic minimum. However, it also 
goes to zero at local minima of {'^ {a)\H\^ [a)) so that 
some caution must be used in interpreting it. 

Note that, unlike with real-time evolution, any errors 
made in integrating the imaginary-time flow equations do 
not accumulate because an accurate step will always take 
the state closer to the ground state irrespective of pre- 
vious steps. The convergence of the energy expectation 
value is quadratic in rj 

so that the approximate ground state energy can be ob- 
tained, to a given precision, with less effort than the 
ground-state expectation value of a general observable. 

2. Conjugate gradient algorithm for finding ground states 

There are a wide range of unconstrained minimization 
algorithms available that often provide far better conver- 
gence than simply taking finite steps along the gradient, 
including the non-linear conjugate-gradient (CG) method 
for approximately quadratic functions (see appendix A). 
Applying such techniques to quantum states restricted to 
a variational manifold may allow us to find ground states 
more efficiently than by integrating the imaginary-time 
TDVP fiow equations. Here we present a naive varia- 
tional implementation of the non-linear CG method that 
can be implemented using only the tools already needed 
for the TDVP. Together with gauge-fixing conditions, it 
is well-defined for uMPS. For more information about the 
differential-geometric properties of A^uMPS, see [29]. For 
more details about optimization on Riemannian mani- 
folds, see [30]. 

The function to minimize is H(x, x) — 
{'i! {x)\HY^ [x]) , which is approximately quadratic 



in the variational parameters x near any stationary 
points. The key difference to the standard CG method 
is the introduction of the non-trivial parameter metric 
gij{x) = {di'^{x)\dj'^{x)). For each step n of the 
algorithm, we require the gradient with respect to a;„, 
which is given by r^^ — g^^ (9i\E'(a;„)|if |5'(a;„)) and which 
we can calculate by minimizing || |<I>(r„)) + H |^(a;„)) ||, 
as with the TDVP. We also need the factor 

Pn — — i 1 ' 

rn-Vn r\^rig,j 

which we can again calculate using methods already 
needed for the TDVP. 

Additional work is required, however, because each it- 
eration of the algorithm involves making a step of length 
a that minimizes H along a given direction Pi . To do this 
in curved space, we should follow a geodesic. Also, to ob- 
tain Pi we must add tangent vectors and be- 
longing to tangent planes at different points Xi and aj^^i, 
requiring the parallel transport of This adds signif- 

icantly to the complexity of the algorithm. However, if 
g{x) is well-behaved such that the parallel- transport map 
is approximately trivial then we can make steps using 
Xn+i sa Xn + anPn with p„ « r„ /?„_ip„_i. Whether 
this assumption is reasonable depends on the particular 
combination of system and variational class. Neverthe- 
less, should it not hold, the line-search used to find the 
step-size still guarantees that the energy will fall with 
each step, so that failure is not catastrophic and merely 
leads to slower convergence. 

In this work, we observe that the above naive method 
is highly effective in the case of uMPS applied to lattice 
(/>'*-theory near its critical point, as exemplified in Fig- 
ure 4. To further improve efficiency, we also implement 
some additional optimizations: For near-critical systems, 
the slowest part of the algorithm, which is also the bottle- 
neck for the TDVP algorithm, is the determination of the 
eigenvectors I and r of E, which we do iteratively. When 
taking small (imaginary) time steps, convergence speed 
improves when using / and r from a nearby state (e.g. the 
previous step) as a starting point for the iteration. In the 
CG algorithm, each evaluation of H for some a visited 
during the line-search requires I and r to be determined. 
To speed this up we store I and r for each point visited, 
using the closest (in terms of a) stored copies as starting 
points for the iteration at each new point visited. Also, 
we do not demand the optimal value of a to high preci- 
sion, since conjugacy will eventually be lost anyway due 
to the assumptions made and because the target function 
is not exactly quadratic. This usually reduces the num- 
ber of evaluations of i? (A, A) to less than ten for each CG 
iteration. Wc use the same optimized line-search routine 
to determine the step size for the gradient-descent results 
in Figure 4. 

We also observe improved convergence of the CG 
method when performing a small number of TDVP steps 
(of fixed step-size) after each reset of the CG algorithm. 
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Figure 4. Convergence of the field expectation value {</>) with 
CPU time for the conjugate gradient (CG) method versus 
imaginary-time evolution via Euler integration of the TDVP 
flow equations and gradient descent (GD — stepping along 
the gradient as with the TDVP, but using a line-search to de- 
termine the size of each step by minimizing the energy). The 
model is (/)*-theory (as defined in section I C) with parameters 
A = 0.2 and X/jl^ ~ 69. The bond-dimension is D = 64 and 
the stopping criterion is 77 < 10~®. The line ("CG final") in- 
dicates the final value taken from the CG curve. We use the 
same line-search algorithm for both the CG and GD meth- 
ods. The gaps in the GD curve are large jumps that could 
occasionally be made in a particular direction. 



C. Excitations with uniform matrix product states 

Given a set of trial states |$(b)) linear in their 
parameters and orthogonal to the ground state, the 
stationary points of the energy functional H{b, b) = 
(<i>(b)|i7|$(6)) / ($(6)|<i>(b)) represent approximate ex- 
cited states. These can be found by solving the gen- 
eralized eigenvalue equation 



H6 = ENb, 



(16) 



where b'TL^th'^ = ($(b)|iJ|$(b')) and h'^stb'^ = 
((f>(6)|$(b')). Given a uMPS approximate ground state, 
the uMPS boosted tangent plane (11) represents a good 
set of ansatz states for probing low-lying excitations of 
uniform systems [15] using this method. 

The suitability of the tangent vectors as ansatz-states 
is based on the ideas of Bijl, Feynman and Cohen and 
assumes that elementary excitations are momentum su- 
perpositions of local disturbances of the ground state. 
Where there is more than one ground state, such as in 
the case of spontaneous symmetry-breaking, elementary 
excitations may also involve their combination to form 
topologically non-trivial states (for example, kink solu- 
tions). For this reason, we additionally include the case 
where a local disturbance interpolates between two de- 



generate ground states. We write the resulting states as 
\^p{B-A,A)) = 



,t 



{4=0 



n 



+00 



n 



(17) 

VR \s) , 



where A = A recovers the boosted tangent vectors for 
uMPS (11) and setting A and A to be the uMPS param- 
eters for two different ground states gives us topologically 
non-trivial excitations. With these ansatz states, which 
are linear in the parameters i?'', excitation energies can 
be obtained by solving (16), which in this case becomes 

UpBi — AEiNpBi, 

where the index i denotes the ith solution, Bi is a vector 
of length dD^ containing the entries of each B^ and the 
matrices Hp and Np are defined as 

2it5{p' - pjB^lipB' = {^p,{B)\H - (H) \%{B')) and 
27r<5(p' -p)StNpB' = {%,{B)\%{B')) , 

where we subtract the ground-state energy (H) = 
{'${A)\H\'^{A)) so as to obtain a finite eigenvalue AEi, 
which is thus the energy difference between the excited 
state and the ground state. 

The effective Hamiltonian term B^HpB' contains three 
(infinite) sums over the lattice sites: One from each 
|$p(i3)) and one from the Hamiltonian. Terms where 
B and B' occur at different lattice sites n and n' ^ n 
acquire a factor e'P("-"'). Infinite sums occur over pow- 
ers of E^, E^, e+'Pi?^ and e'^^E^, where the first 
two have spectral radius 1 and can be calculated with 
techniques used in the TDVP algorithm (section II B), 
leading to pseudo-inverse factors {I - E^f. Ef and 

E^ are related to the overlap between the two ground 
states. The per-site fidelity is equal to the spectral ra- 
dius p(i?^) = p{E^) which, unless the states are the 
same (up to a phase), is less than one. For two differing 
ground states, these infinite sums thus become geometric 
series Y.n=oi^^'^ ^A^"" = (I " e+'PEj)-\ If the states 
are the same, then p{E^) = 1 and the inverses must 
be replaced by pseudo-inverses. For example, a part of 
B^HpB' wh 
separated is 



B^HpB' where all three summed-over lattice sites are 



+00 +00 



J2 E miiEjr-'Ef iEjr'''Hii\f) 



m— 1 m' — l 



{1\E^{1 - e+'PEjr'Efil EifHii\f) , 



where r is the right eigenvector of E^, m is the number 
of sites between B and B' , m! is the number of sites 
between B' and the Hamiltonian term h and we assume 
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Hit 

AA 



< 1. The Hamiltonian term is contained within 



ttAB 
J^CD 



stuv 



st\h ~ (h) \uv) A'B* ® C'D". 



There is no infinite term corresponding to the pseudo- 
inverse in the above example because (llH'^^lf) = 0. 
Additional simplifications can be made by implementing 
the gauge-fixing condition 



If) 



= 0. 



factor) or with e"'"'^ (e ^p) (resulting in the momentum 
shift). We adhere to the convention of [15] and demand 
that the largest eigenvalue of is real and positive 
which, in the case of equivalent states differing only by a 
phase A = e'*yl, corresponds to = 0. 



1. Mean-field case 

In the mean field case D = 1, where there is no inter- 
site entanglement, the uMPS excitation ansatz simplifies 
further. The trial states are 



A corresponding parametrization of is 



(18) 



where x G -A^DxD(ti-i) £^nd the D{d — 1) x dD matrix 
[y]{a,s);p = IS defined so that contains an or- 

thonormal basis {VV^ = I) for the null-space of R\ with 

resulting in VR = 0. For A = A, this parametrization 
is identical to (13). With it, the overlap term becomes 
B{xY'NpB{y) = tr[a;'''y] = {x\y) so that the problem 
turns into a standard eigenvalue problem. The effective 
Hamiltonian term becomes 

B\x)n,B{y) = m^lUp) + 



{l\E 



B{x)A 



'AB{x) I 



B{x) 



"^ABlx 



r) 

AAi^ 



B{x)\ 

^aa"^) 

+ e~-^^ {l\E"^^^\l e--^^Ei)-'Ei^^-^{l E^f Hll\f) 



e'- {l\E^^^^{l-e+'PE^)-^Efy\l 



EifHj^^ 



s+'^(ZK,)(I-e+'^'£;j)-ii/|J)^|r~) 



f) 



-'P {l\El'^y\l-e-'PEi)-^H^^ 



B(x)A 



+ e+^'P {l\E^^^^{l-e+'PE\)-'H^''^y\f) 



{l\Ef'"\l - e-'PEi)-^H 



AA 

AA 
AB(x) 



(19) 



where we again note that the inverses turn to pseudo- 
inverses if ^ = A. It is possible to implement these 
operations with 0{D^) time complexity, avoiding direct 
calculation of inverses as in the TDVP algorithm (see 
section II B). A sparse eigenvalue solver can then be used 
to efficiently obtain eigenvalues. 

A final ingredient is needed to define the momentum 
p in the case A ^ A, because an overall phase on A 
effectively shifts the momentum of the ansatz states 

\<^p{B-e'^A,e'^A)) ^ \%+^^^{B- A, A)) , 

which can be seen in (19), where every factor A (A) is 
paired either with (A^) (cancelling any extra phase 



\%{b:a,a)) ^ 



IV'(a)) 



(E> \ipib)) (g) 



+ 00 



|V(a)) 



.n+l 



where the rank 3 tensors A and B of (17) have become 
vectors a,b G C^. In this case, the x operators E 
and the corresponding vectors are just numbers so that 
the requirement p{E^) — 1 with the largest eigenvalue 
being 1 implies E^ = a • a = 1. The normalized "eigen- 
vectors" {l\ and |r) are thus also equal to 1 and projecting 
them out using Q = 1 — |r) (^| leaves zero. All terms in 

(19) containing the pseudo-inverse oil - Ei ov I ~ Ej^ 

thus drop out. In the case A — A, this leaves only the first 
four terms. Otherwise the inverse factors il-e+'PEj^)-^ 

and (I — e ^^^E^)^^ are just positive numbers and the 
last four terms are non-zero as well. 

The gauge-fixing condition B'^A'^'^ = corresponds 
to the elimination of the norm degree of freedom, where 
the parametrization (18) forces b into the subspace or- 
thogonal to a. We can directly obtain the effective 
Hamiltonian as a d x d matrix 



[H 



s^\h'\'ipt) +e~'P {'il;s\h'\ti^) 



+ e-'P (V-li) (1 - e-'P (Vl^))"' {s^\h'\^i>) 

(^it)(i-e-'p {mr^ {Mh^i'i^)] , 



+ {s\i>) (1 



+ e 



-2ip 



where = |V'(^)) a-nd h' = h— (h). In the topologically 
trivial case, where a = a, the terms in square brackets 
drop out due to b ■ a — b ■ a — 0. 



III. STUDYING QUANTUM FIELDS WITH 
MATRIX PRODUCT STATES 

In this section, we use the variational conjugate- 
gradient method for uniform matrix product states 
(uMPS) of section II B 2 to determine the continuum 
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critical parameter of (/)^-theory, improving on previous 
numerical results. We also study the special case of uni- 
form mean-field theory (MFT) states, which correspond 
to uMPS with bond dimension one. As well as the vac- 
uum expectation value of the field, which plays the role of 
the order-parameter (see section ID), we investigate the 
energy of the lowest-lying excitation as a phase-change 
indicator, which tends to zero at the critical point and, 
by universal correspondence to the Ising model, should 
scale linearly in its vicinity. Furthermore, we extract the 
central charge of the conformal field theory (CFT) of the 
critical system [31], which is also expected to be univer- 
sal [32-34] and calculate the spectral density function of 
the near-critical lattice theory. 



A. Method 

As set out in section IC, (1 + 1) dimensional 0^-theory 
can be put on a spatial lattice, in a way that formally 
recovers the continuum theory in the limit of zero lattice 
spacing a — >■ 0, using the nearest-neighbour Hamiltonian 



A 

4!' 



^^a? are dimensionless pa- 



where A = Xa? and 
rameters. The theory exhibits spontaneous symmetry- 
breaking, as detailed in section ID, where a particular 
value of A//2|j(A) characterizes the critical point for a par- 
ticular lattice-spacing a, hence the dependency on A(a). 
/i/f = //q -|- 5jl\ is the renormalized mass, which is finite 
for a > and is given by (5). We use the uMPS conju- 
gate gradient algorithm of section II B 2 to obtain ground 
states up to some tolerance 77 (see section II Bl) giving 
us access to approximate ground-state expectation val- 
ues, and the uMPS excitation ansatz of section II C to 
obtain excitation energies. To study the system using 
uMPS, we first need an appropriate basis. 



1. Position basis with a cut-off 

To represent states using the uMPS formalism we 
choose the position basis described in section IC: 



|0„) [a„, aJ,J = 6n 



(20) 



1 
71 



V2 



(4 



We provide the matrix-elements of relevant operators for 
this basis in appendix B. Since the site subspace is infi- 
nite, we must introduce a cut-off so that states can be 
stored using a finite number of parameters. We there- 
fore limit ourselves to Hn = C'' such that the highest 
available number-eigenstate is \d— 1), assuming that a 
good approximation to the ground state does not require 
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Figure 5. Illustration of the field shift needed to centre fluc- 
tuations about zero. 



the higher modes to be present. That this should be the 
case for the symmetric phase seems intuitive consider- 
ing the form of the classical effective potential (see Fig- 
ure 2), but things are less clear for the symmetry-broken 
case where the ground state is centred about one of two 
separated wells away from the origin. The cut-off may 
thus affect the accuracy of symmetry-broken states more 
significantly than symmetric ones. We also expect the 
higher modes to be more important for states near to 
the critical point, where fluctuations diverge. 



2. Field- shifted basis 

It should be possible to avoid higher excitations in the 
symmetry-broken phase [cj)) 7^ 0, thus mitigating the ef- 
fects of the cut-off, by changing the basis such that the 
operator 0' in the new basis has an expectation value of 
approximately zero {(j)') « 0. We effectively shift the ori- 
gin in a plot of the effective potential by some amount 
towards the minimum, such that fluctuations are cen- 
tred about (j)' = 0. That higher excitations in the shifted 
number basis are then avoided seems intuitively reason- 
able given the classical effective potential, where each of 
the two wells (in the symmetry-broken case) looks locally 
similar to a single- well potential. Figure 5 illustrates this 
procedure. 

The change of basis corresponds to the unitary 

f/(0c) = e''^=^ 

with TT being the conjugate momentum operator of (20). 
It defines new creation and annihilation operators 



such that <p'^ = 4>n ~ (f'c where G K characterizes 
the shift. In terms of operators in the shifted basis, the 
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Figure 6. Histogram plots for tlie number operator in the 
shifted basis at varying distances A//i|j = 67,100,200 from 
the critical point X/p-^,^ ~ 66 in the symmetry-broken phase 
(with A = 0.1 and D = 128,64,64 respectively). The higher 
modes carry more weight for states nearer the critical point. 



Hamiltonian is 



H 



E 



Using this Hamiltonian with a value of (j)c ~ (^) should 
thus help to avoid higher excitations and allow us to effi- 
ciently represent ground states with large values of | (0) | . 
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Figure 7. Histogram plots for the number operator in the 
shifted and non-shifted bases near (top, X/ji\ = 80) and 
far from (bottom, A//i|j — 200) the critical point in the 
symmetry-broken phase (A = 0.1). The effect of shifting by 
approximately {(f)) is much stronger far into the symmetry- 
broken region (higher (0)), where we see a weight-shift from 
the higher modes to the zero mode. The shifted states were 
obtained with d — 16, the non-shifted with d = 24, hence the 
greater range of the non-shifted points. All four states have 
D = 64. 



3. Effects of the Hilbert space cut-off 

The effects of the local Hilbert space cut-off d are, as 
expected, relatively strong near to the critical point, be- 
coming weaker further into the symmetry-broken phase 
when using the shifted basis (see Figure 6). Without 
the basis shift, states with large values of (0) exhibit a 
weight-shift towards higher modes, as illustrated in Fig- 
ure 7. In all cases, excitation of higher modes drops off 
exponentially, with d = 16 being sufficient to capture the 
most significant contributions, as demonstrated in Figure 
8. 

The shifted basis has an added benefit when sweep- 
ing in the broken phase and using the previous 
ground state approximation as a starting state for the 
next ground-state search. In this case, adjusting the shift 

towards the next predicted ((/>) (according to a prelim- 
inary fit of (21)) improves the starting state by bringing 
{(j)) closer to the new ground-state value, leading to faster 
convergence. This is because a shift of (0) always centres 
the state about the origin in the shifted basis (see Figure 
9). Adjusting the shift by some also adjusts {(j)) by 
the same amount. 



4.. Locating the critical point using the field expectation 
value 



As noted in section ID, since (t^) is the order- 
parameter associated with the (/("^-theory phase-change, 
it can be used to identify the critical point. A possible 
strategy for finding the critical parameters for a, A > 
might thus be to fix A and sweep until one sees a tran- 
sition from (0) 7^ to {(j)) ~ or vice versa. However, 
this is not practical because the amount of entanglement 
in the ground state (for example, as quantified by the 
half-chain entropy (9)) tends to infinity as the critical 
point is approached, such that accurate representation 
using uMPS requires the bond-dimension D to approach 
infinity also. Since the computational complexity of the 
TDVP algorithm scales as 0{D^), this bisection method 
cannot achieve high accuracy for reasons of practicality. 

Instead, we approach the critical point from the 
symmetry-broken phase, noting that physical quantities 
obey power laws in the vicinity of critical points (see sec- 
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Figure 8. Scaling of (0) (blue) and the half-chain entropy S 
(red) with the Hilbert space cut-off d far into the symmetry- 
broken phase (top, = 200, D = 64) and near to the 
critical point (bottom, = 67, D — 128) using a shifted 
basis. In both cases A — 0.1. It appears that d = 16 is 
sufficient both near and far from the critical point. Additional 
variation for d > 16 in the near-critical case is due to high 
sensitivity to the level of convergence (states were obtained 
with a tolerance of < 3 ■ 10~^). 



tion ID). For (0) we can thus write 

A 



A 



(21) 



where A{X) is a constant, /3(A) is the critical exponent 
and flj^ ^(A) is the critical value of /ifj for a given A. Fit- 
ting this equation to (0) as a function of A/ p,^ (with fixed 
A) as near as possible to the phase-transition, we obtain 
an estimate for the lattice critical parameter A//2|jj,(A). 

We can then use a series of fits with A ^ to extrapo- 
late an estimate for the critical parameter A//i^ of the 
continuum theory. 

Initial simulations show that, as expected, the half- 
chain entropy S of the ground state approximation tends 
to infinity as the critical point is approached. This is vis- 
ible in Figure 10, where we show results obtained from 
high bond dimension limits as well as using a fixed bond- 
dimension. As further confirmed in Figures 11 and 12, 
a fixed D is not sufficient to capture ground states near 
the critical point. Note that a phase transition does. 
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A/Afl 




Figure 9. Visualization of the basis shift far into the 
symmetry-broken phase (A = 0.1, A//if{ = 200). histograms 
are plotted with (green) and without (blue) the basis shift 
(j)^ ^ {(f>) . Additionally, the histogram of the shifted-basis op- 
erator (ji' is plotted for the shifted state (red). States were 
obtained with d — 16, D = 64. 



Figure 10. An example plot of the order parameter (0) 
(circles) for fixed A — 0.5, sweeping \/fL\. The half-chain 
entropy 5* is also shown (triangles). The cyan and red curves 
represent high bond-dimension limits with D < 80, whereas 
the blue and magenta curves are for fixed D = 32. All ground 
state approximations are converged to a state tolerance rj < 
lO"''. 

in fact, occur for fixed D, albeit not at the exact crit- 
ical point, but at increasingly lower values of X/jlj^ for 
decreasing values of D (and for decreasing A). This is 
consistent with the entropy shown in Figure 10, which 
is asymmetric about the critical point, falling off more 
slowly in the symmetric phase. Since a fixed D repre- 
sents an upper bound on the amount of entanglement in 
the state (see section II A), the uMPS variational mani- 
fold 7MuMPS comes closer to the exact ground-state when 
its entropy S is lower. Given the asymmetric entropy 
of 0^-theory, this implies that symmetry-broken ground- 
states are easier to approximate than symmetrical ones 
(for a given distance in parameter space from the critical 
point). For a symmetric ground state with high entropy. 
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Figure 11. A plot of (0) for fixed A — 0.5, sweeping 
for several bond dimensions. At lower values of D, finite- 
entanglement effects shift the apparent location of the critical 
point to lower values of \/fl\. 




Figure 13. Illustration of the relationship between the uMPS 
variational manifold for fixed bond-dimension AIumps and the 
ground state \Q,) for parameters close to the critical point in 
the symmetric phase. Hilbert space is divided into symmet- 
rical and asymmetrical (in (j)) states. 
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Figure 12. Scaling of the half-chain entropy S (top) and of 
{(p) (bottom) with the logarithm of the bond dimension D for 
points near (green, — 70) and far (blue, A//i|j — 196) 

from the critical point (A = 0.1). A higher bond dimension 
is necessary to accurately represent near-critical states com- 
pared to far-from-critical states. 



a low-lying excited state with much lower entropy may 
thus turn out to be the best available ground state ap- 
proximation in A^uMPS- Such an excitation should be 
available for such states, since a small change in the pa- 
rameter jj,^ results in an asymmetric ground state (on 
the other side of the critical point). The situation is il- 
lustrated in Figure 13. The same Z?-dependent shift of 



the critical point is observed with the transverse Ising 
model [35]. 

One might consider using this behaviour together with 
finite-entanglcmcnt scaling techniques [35] to obtain in- 
formation about the true critical point (for example, the 
critical exponent), but this requires precise knowledge of 
its location. Instead, we take data at several values of D 
in order to obtain high-£) limits of (i^), which we then fit 
using (21) to obtain an estimate for the location as well 
as the critical exponent. 

Since we take a high-D limit of the approximate 
ground-state value of {(p) for each parameter combina- 
tion (requiring a higher D for the higher-entropy states 
closer to the critical point) there is a practical limit on 
how near we can come. This is unfortunate, since the 
fit (21) is highly sensitive to near-critical points, where 
the gradient goes to infinity. Power-law scaling is also 
only exactly fulfilled infinitesimally close to the critical 
point, such that including data points further away de- 
creases accuracy. We thus only fit the points closest to 
the critical point for which we have sufficient (in terms 
of bond-dimension) data. 



5. Locating the critical point using excitations 

Another approach to finding the critical parameters, 
given A, is to plot the energy of the lowest-lying excitation 
AE = aAE against X/jij^, which should tend to zero 
as we approach the critical point X/fifj cW from either 
side. We describe an ansatz for obtaining the lowest-lying 
excitation energies, given a uMPS approximation to the 
ground state, in section II C. 

In obtaining excitation energies in the symmetry- 
broken phase, topologically non-trivial excitations must 
be taken into account. This precludes the use of the 
(/)-shifted basis mentioned in section III A 1 to represent 
the state, since approximations to both possible ground 
states are required for the calculations and these must use 
the same basis. For reasons of efficiency, it thus makes 
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sense to focus on states near to the critical point where 
the shifted basis is not needed. This should not cause 
problems since this is where we expect power-law scaling 
to be more exactly fulfilled. 

We first locate the lowest-lying excitation of the 
symmetry-broken phase and determine whether it is 
topologically trivial or non-trivial, whilst confirming that 
it goes to zero for some value of A//i^. To do this, we 
use the excitation ansatz to determine dispersion rela- 
tions for the lowest-lying topologically trivial and non- 
trivial excitations for fixed A and several values of X/jj.%- 
We then use linear extrapolation of the excitation ener- 
gies at each momentum to obtain a dispersion relation 
at the first point where one of them goes to zero, which 
should correspond to the lattice critical point A//i|j ^.(A). 
The result is shown in Figure 14, where we see that the 
lowest-lying excitations are the topologically non-trivial 
soliton (kink) excitations (at zero momentum). A plot 
of this excitation energy versus A//i^ exhibits almost ex- 
actly linear scaling, consistent with the transverse Ising 
model, suggesting the use of linear regression to obtain 
an estimate for the critical parameter. The plot is shown 
in Figure 15, which contains the excitation energies ob- 
tained for several bond-dimensions. The small change in 
excitation energy near the critical point when increas- 
ing the bond-dimension from D — 16 to D = 48, com- 
pared with the change in (0) shown in Figure 11 (for a 
larger lattice-spacing), suggests that finite-entanglement 
effects are less severe for the excitation energy than for 
{(j)). Certainly, the ground-state energy should reach a 
high-D limit sooner than (0) simply because the approx- 
imate ground state is close to the energy minimum. Also, 
if the exact lowest-lying excitation is highly localized, we 
should need only a relatively low D to approximate it 
well. It seems excitations present a less computationally- 
intensive way of obtaining a good estimate for the critical 
parameter, compared with (0). 

We note that, for fixed D, the lowest-lying (soliton) 
excitation receives a negative energy (with respect to the 
approximate uniform ground-state) for sufficiently low 
values of A//i|. (see Figure 15). The position of this 
crossover is very close to the critical point predicted by 
linear extrapolation. The existence of negative approxi- 
mate excitation energies indicates that the topologically 
non-trivial ansatz states include a better approximation 
to the exact ground state than A^uMPS- This is consis- 
tent with the Z3-dependent shift of the apparent phase- 
transition of ((/)): If the exact ground state is symmetric, 
but the approximate ground state is asymmetric, a topo- 
logically non-trivial excitation interpolating between the 
two degenerate asymmetric approximate ground states 
should be locally closer to the exact ground state at the 
disturbance. 

In fact, we can construct a better uMPS ground-state 
approximation |^'(A')) using negative-energy kink "exci- 
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Figure 14. Extrapolated dispersion relation at tlie approx- 
imate lattice critical point — 1.0) ~ 64.4 showing 
the lowest-lying topologically trivial (TT) and topologically 
non-trivial (NTT) excitations. The zoomed area shows that 
the non-trivial excitation is the lowest-lying excitation at zero 
momentum. Momenta Q < p < n / a are shown using p = ap. 
Energies l\E are relative to the approximate ground-state en- 
ergy. The bond-dimension is D = 32 and points were linearly 
extrapolated from data at A//i|{ — 70, 75, 80. 
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Figure 15. A parameter sweep of the energy of the lowest- 
lying excitation, which is a soliton (with zero momentum) for 
A = 0.2 at D = 16, 32, 48. That the energy becomes negative 
for low values of indicates that these points lie in the 

symmetric phase (see main text). The line represents a fit to 
the data with _D = 48 using points X/ij-% = 67. . . 70. The 
fitted value of the critical parameter is A//i|j ^ — 65.82(1). 



tations" by defining new 2D x 2D parameter matrices 



eW A' 



where A and A are the parameters for the two original 
ground states, B and B are the tangent-vector parame- 
ters for a kink and an anti-kink and e e M. The resulting 
state contains the original ground states as well as kink 
states at order e plus multi-kink states at O(e^). This 
leads to energy contributions e^i?kink (there are no kink 
contributions at order e because the kinks are orthogonal 
to the original ground states). Since -Ekink < ^i'i'(yi)) the 
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state l^'(A')) can, depending on e, have a lower energy 
than e can be interpreted as a kink-density with 

an optimal value depending on the higher-order energy 
contributions. 



6. Mean-field theory 

When D = 1, the uMPS variational class is the same as 
that of the uniform product states (or mean field theory 
(MFT) states) 

|*(a)) = • • • (g) |?A(a)) (g) IV'(a)) (g . . . 



where \^{a)) = J2s=o 1^) ^^"^ ^ ^ ^^^^ case, the 

per-site energy expectation value takes on an effective 
one-particle form 



2 2 ^ 4!^ 



= (V'(a)l 



\^Pia))+alia), 



where ct^ is the 0-variance a^{a) — (i/)(a)|(/)^|?/;(a)) — 

{ip{a)\(t>\ipia))'^ . An approximation to the ground state 
can then be found by applying the time-independent vari- 
ational principle and minimizing (h) with respect to the 
d parameters a, which can be taken to be real since all 
matrix elements in the above expression are real. The 
gradient of {h) is also readily obtainable 
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making many commonly-used minimizing algorithms ap- 
plicable, such as the quasi-Newton method of Broy- 
den, Fletcher, Goldfarb, and Shanno (BFGS) [36]. This 
method is much simpler and more efficient than apply- 
ing the imaginary-time TDVP algorithm for uMPS with 
D = l. 

Normalization presents a minor complication. The 
above equations assume {ip(a)\'ip{a)) = 1, which im- 
poses a constraint = 1 on the variational parameters. 
Rather than using a constrained optimizer, we eliminate 
the norm degree of freedom by switching to d-dimensional 
spherical coordinates such that the norm corresponds to 
a single parameter and can easily be fixed and ignored. 

To estimate the location of the lattice critical point 
A/M/jc using MFT, we again obtain ground states for a 

sweep of A//i|j for some fixed A. As in the more general 
case of uMPS with fixed D, we are putting a restriction 
on entanglement by using MFT (D = 1) and thus expect 
an apparent phase-transition to occur at some value of 
^/P-R < ^/P-'r c fo'^ ^ given A. We use the location of the 
apparent transition as an estimate for X/flj^^, obtaining 



it from both {(j>) and excitation energies calculated us- 
ing the MFT excitation ansatz of section II CI. Since 
we do not expect power-law scaling of physical quanti- 
ties to be reproduced by MFT, we do not attempt to 
fit data using power laws. Instead, we use bisection to 
pin down the apparent phase-transition in (j), which is 
possible due to the relative ease of finding MFT ground 
states, and interpolate the lowest-lying excitation ener- 
gies (in the apparently symmetry-broken phase) to obtain 
the point at which they become negative (these excita- 
tions are topologically non-trivial, as for D > 1 — see 
the above explanation). 

Although we can obtain estimates for the lattice criti- 
cal point using these methods, we do not expect to obtain 
useful information about the continuum critical theory 
due to the lack of entanglement. However, since it is also 
possible to interpret the lattice critical point as a con- 
tinuum limit of a non-critical theory (see section ID), 
an ability to estimate its location using mean-field the- 
ory indicates that useful predictions about non-critical 
continuum theories can be made. 



7. Central charge 

We determine the central charge associated with the 
conformal field theory of the critical system using finite- 
entanglement scaling techniques. It is known that, for 
infinite one-dimensional systems with a second-order 
phase-transition, the half-chain entropy of the ground 
state in the vicinity of a critical point with conformal 
invariance is 



S 



6 



log(e/«), 



(22) 



where ^ is the correlation length and c is the "central 
charge" [33]. Approaching the critical point, ^ — oo 
and the entropy diverges. The central charge specifies a 
conformal field theory (CFT), which describes behaviour 
at the critical point in the continuum limit. 

We know from (9) that the maximum half-chain en- 
tropy of a uMPS state, contained, assuming right canon- 
ical form, in the diagonal entries of the D x D matrix 
I, is directly related to the bond-dimension D, which 
is also the maximum Schmidt-rank of the correspond- 
ing Schmidt decomposition. Thus, for lower values of 
D, finite-entanglement effects occur and the value of S 
scales with D. It turns out there is a simple relationship 
between S, D, and c describing this scaling [31] 



S = 



1 



yi2A 



1 



log I?. 



(23) 



We can thus obtain an estimate for c from values of 
S taken from a number of ground state approxima- 
tions with varying D. For this to work, we must be 
close enough to the critical point so that (22) is valid 
and use small enough D so that S is limited by finite- 
entanglement effects. We can then use linear regression 
to fit (23) and obtain c. 
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B. Results and analysis 

1. Estimates of the continuum critical parameter 

Figure 16 shows estimates for the critical parameter 
2 (A) taken from sweep plots of (0) and of the lowest- 
lying excitation energy Ai?, approaching the continuum 
limit A — )• 0. 
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Figure 16. Approximate values for the lattice critical param- 
eter A//i|j^(A) (top) and the (0) critical exponent /3(A) (bot- 
tom) obtained from linear fits to the lowest-lying excitation 
energies and from power-law fits to the order-parameter 
{4>) for values of A approaching the continuum limit A — >■ 0. 
The line in (a) corresponds to the fourth fit of Table I. 

The two sets of values show good agreement, with the 
largest discrepancy occurring for A = 6, where we found 
high-Z? limits of (0) particularly close to the critical point 
without resorting to very high bond-dimensions. Exclud- 
ing the points of lowest (</>) from the fit pushes the fitted 
value of X/fijic upwards, closer to the AE value, lead- 
ing us to speculate that the excluded (</<) values were not 
accurate enough, possibly due to insufficient convergence 
of the uMPS ground state. We are inclined to trust the 
results obtained from the AE data over those from fits 
to ((/)), in particular due to the relative robustness of the 
linear regression fit to errors made near the critical point. 

As expected (see section ID), non-linear behaviour of 
A//i^^(A) is present. Given that the exact behaviour is 
unknown, but is predicted to be logarithmic, we follow 
[21] and fit a series of functions, evaluating the statis- 
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h ciA 


+ cgAln A H 


- C3 A^ In A 


66.26(3) 


24.8 


66.42(5) 


6.22 



Table I. Fits, for lattice-spacings approaching zero, of the 
lattice critical parameter ^(A) obtained from power-law 

fits to uMPS ground-state ((;/))-values and from linear ex- 
trapolation of lowest-level excitation energies AE (all in the 
symmetry-broken phase), /c = X/nji^c is the extrapolated 
continuum critical parameter. We limit the {(j)} data fitted to 
obtain each X/jj-'ji^c to a few points close to the critical point 
with {(f>) < 0.59. The fitted data is plotted in Figure 16. 
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Figure 17. Position in \/ fiji of the apparent phase-transition 
obtained from the mean-field results for {(f)) (green) and the 
lowest-lying excitation energy AE (red) compared to the crit- 
ical parameters X/jl\ ^(A) obtained from the uMPS AE data 
(blue). 



tic to judge which can be reasonably used to predict a 
continuum value A/^|j ^. The results of the fits are listed 
in Table I, where we define our final estimates for A//i|j. ^ 
to be the fitted values with reduced statistic x^/dof 
closest to one. 

We find that the critical exponent /3(A) obtained only 
from fits to {(j)) agrees poorly with the predicted trans- 
verse Ising value of 0.125, the fitted values near the con- 
tinuum limit being significantly higher, as shown in Fig- 
ure 16. This we attribute to insufficient data near to the 
lattice critical points, noting that the effect of excluding 
the points of lowest {(()) is to increase the fitted value of 
/3(A) further. Using the A//i|j ^(A) values taken from the 
AE' data together with the (</)) data, we obtain a second 
estimate of /3(A) that, in the continuum limit A -> 0, 
shows a much clearer trend towards the Ising value, in 
support of the greater reliability of the Ai5-based esti- 
mates of the critical parameter. 
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2. Mean-field results 

Due to the lack of entanglement in the mean-field class 
(see section III A 4) , we expect a large shift of the appar- 
ent phase-transition to lower values of A//i|.. We ob- 
serve that the shift takes the apparent phase-transition 
(in ((/)) and in AE) towards A///^ = in the contin- 
uum limit A — > (see Figure 17). For higher values 
of A, the apparent phase-transition starts to approach 
the approximate critical parameters obtained above us- 
ing uMPS with D > 1. The mean-field excitations data 
gives better predictions for the critical parameters than 
{4>), in agreement with the relatively low sensitivity of 
the excitation energies to the bond-dimension observed 
with uMPS. 



3. Central charge 



We extract the central charge by finding approximate 
uMPS ground states for various low bond-dimensions at 
(or very near) the lattice critical point as determined 
from the uMPS excitations data. We observe approxi- 
mately linear scaling of the entropy, as predicted by (23), 
with the central charge extracted from the gradient of a 
linear fit agreeing well, for larger values of A, with the 
prediction of c = 0.5 from the transverse Ising model. 
For lower values of A, we find that the gradient often 
agrees poorly with c — 0.5, despite the scaling remaining 
linear. At this point, we do not have a good explanation 
for this discrepancy, so we leave it as a subject for future 
investigations. Our results are summarized in Figure 18. 



4- Spectral density 



To further demonstrate the convenience of having ap- 
proximate ground states in uMPS form and the useful- 
ness of the uMPS excitation ansatz, we obtain the spec- 
tral density function (3) from the overlap of approximate 
excited states with the approximate uMPS ground state. 
Results for states in the symmetry-broken and in the 
symmetric phase are shown in Figure 19 and Figure 20 
respectively. By interpreting the lattice critical point as 
a continuum limit of a non-critical theory (see section 
ID), a series of such plots could be used to extrapolate 
a continuum spectral density for that theory. The na- 
ture of the uMPS excitation ansatz means that not all 
excitations can be captured. This explains the lack of a 
continuum of excitations, expected soon after the single- 
particle state. It would be interesting to compare the 
results with those of other approaches such as [37]. 



CO 




A c 



0.1 66.2 0.372(15) 
0.2 65.81 0.375(14) 
1 64.8 0.419(13) 
3 64.5 0.495(23) 

5 64.9 0.475(19) 

6 65.1 0.508(22) 



0.109(9) 

0.086(7) 

0.019(76) 

-0.029(10) 

-0.041(9) 

-0.051(10) 



Figure 18. Entropy-scaling with D for the near-critical lat- 
tice theory, using estimates for obtained from ex- 
citations data. The legend shows the value of A. The ta- 
ble shows the parameters used and the values for the critical 
charge c and the S-intercept a derived from the linear fits. 
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Figure 19. Spectral density function in the symmetry-broken 
phase (A = 1, A//l|j = 77) obtained at D = 16. Dirac delta 
functions are replaced by Gaussians to aid visualization. 



C. Discussion 

The consistency of the lattice critical parameters ob- 
tained from approximate uMPS ground-state field expec- 
tation values {(f) and those derived from the lowest-level 
excitation energies Ai? calculated using the uMPS ex- 
citation ansatz, as well as agreement of the critical ex- 
ponents with their transverse Ising counterparts, demon- 
strates the validity of both methods for studying critical 
(1 + l)-dimcnsional (^"^-thcory. The finite-entanglement 
scaling method we use to obtain estimates for the cen- 
tral charge c of critical (/)^-theory shows promise: Ap- 
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Figure 20. Spectral density function in the symmetric phase 
(A — 1, X/jl^ = 60) obtained at D = 48. Dirac delta functions 
are replaced by Gaussians to aid visualization. 



proximate linear scaling of the entropy with log2(I?) is 
observed, as predicted, and the fitted values for c agree 
with the transverse Ising value c — 0.5 for higher val- 
ues of A. However, further work is needed to explain the 
discrepancies observed for lower A. 

We find uMPS to be an excellent class of ansatz states 
for studying the critical phenomena of (l + l)-dimensional 
(/)*-theory due to the amount of entanglement in near- 
critical ground states being the main barrier to their 
efficient representation. With MPS, the amount of en- 
tanglement that can be represented is controllable via 
the bond-dimension D, which can easily be varied to 
obtain the limiting behaviour of quantities such as (0). 
In this way, we can avoid errors originating from finite- 
entanglement effects. Additionally, working directly in 
the thermodynamic limit of an infinite lattice completely 
avoids additional finite-size effects and the need for fur- 
ther scaling investigations. Obtaining ground states us- 
ing our variational conjugate gradient method (or the 
TDVP with imaginary time evolution) is also very con- 
venient, since we need only enough storage capacity to 
capture the approximate state at one point in time, un- 
like when simulating the Euclidean theory on a space- 
time lattice. Also, using the TDVP, the computational 
complexity scales linearly in t. 

The class of uniform mean-field states (uMPS with 
D = 1), despite featuring no entanglement, appears 
suited to estimating properties of non-critical continuum 
theories in some cases, even if not of critical theories. 
Owing to the relative ease with which MET ground-state 
approximations and excitation energies can be obtained, 
and the low computational cost of extending them to 
higher space-time dimensions (due to the lack of entan- 
glement), these methods represent another useful tool for 
investigating quantum field theories. 



1. Comparison to other methods 



We summarize existing literature estimates for the con- 
tinuum critical parameter A//i|j^ in Table II, where we 
include our results from Table I with x^/dof values clos- 
est to one. Our estimates agree poorly with the DMRG 
result of [18], but relatively well with the Monte Carlo 
results of [21]. With regard to the DMRG results, where 
the technique used is similar to ours, we can attribute 
the large difference to finite-entanglement effects, since 
these serve to shift the apparent critical point to lower 
values of with a larger shift for smaller lattice- 

spacings. The DMRG parameters used in [18] correspond 
to D — d — 10 [38], resulting in a relatively large shift 
(see Figure 11). The DMRG study also uses only two lat- 
tice critical points to extrapolate a continuum value and, 
as such, misses the non-linear behaviour of A//2|j^(A). 
The Monte Carlo methods used in [21] are very differ- 
ent to ours. They work with the Euclidean theory on 
a finite two-dimensional lattice whereas we work on an 
infinite spatial lattice in continuous time (numerical in- 
tegration of the TDVP fiow equations could be seen as 
analogous to discretizing imaginary time on a lattice, in 
which case our temporal "lattice" is of the length neces- 
sary to produce sufficient convergence of the approximate 
ground-state). Rather than taking finite-size scaling lim- 
its, we take finite-entanglement scaling limits to obtain 
our ground-state approximations. Noting these differ- 
ences, the fact that our results agree to within 2% gives 
us confidence in the methods used. 



Method 


fc 


Reference 


uMPS, TDVP, AE 


66.46(5) 


This work 


uMPS, TDVP, {(P) 


66.30(2) 


This work 


Monte Carlo 


64.8«:«3 


[21] 


Gaussian effective potential 


61.632 


[39] 


Gaussian effective potential 


61.266 


[16] 


GEP and oscillator rep. 


61.26 


[40] 


Spherical field theory 


60.3 


[41] 


Diffusion Monte Carlo 


60 ±4.8 ±2.4 


[42] 


DMRG 


59.89(1) 


[18] 


Continuum light-front 


59.46 


[43] 


Connected Green function 


58.70 


[39] 


Coupled cluster expansion 


22.8 < /c < 51.6 


[44] 


Discretized light-front 


43.95, 46.26 


[45] 


Discretized light-front 


43.70, 33.00 


[46, 47] 


Random phase approximation 43.2 


[48] 


Non-Gaussian variational 


41.28 


[49] 



Table II. Summary of results for the continuum critical pa- 
rameter fc = X/u'r^c from the literature, including our results 
derived from lowest-lying excitation energies AE and from 
{(p), where we use the results corresponding to the x^/dof 
values closest to one (see Table I). 
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CONCLUSION 

The class of uniform matrix product states (uMPS) ap- 
pears well-suited to the study of critical quantum fields 
in (1-1-1) dimensions via lattice regularization. Using 
variational methods like our naive variational conjugate 
gradient method or the imaginary-time time-dependent 
variational principle (TDVP) (where the former provides 
significantly improved convergence speed for the system 
studied), good approximations to ground states can be 
obtained efficiently, even near to the critical point. Here, 
the correspondence between the bond dimension D and 
the maximum entanglement of a state allows the use of 
finite-entanglement scaling to judge the accuracy of phys- 



ical quantities calculated. Compared to Monte Carlo sim- 
ulations, uMPS allows us to work directly in the thermo- 
dynamic limit and has storage requirements independent 
of the imaginary time dimension. Further, low-lying exci- 
tation energies are straightforward to calculate, enabling 
the study of dispersion relations and the spectral density. 
Even mean field theory shows potential for delivering use- 
ful predictions about non-critical continuum theories. 
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Appendix A: Minimization using conjugate-gradient 
methods 



When minimizing func- 
tions that are approximately 
quadratic in their parameters, 
which is always true in the 
vicinity of a minimum (assum- 
ing sufficient differentiability), 
making steps along the gradi- 
ent direction (gradient- descent) 
is often a sub-optimal way of 
reaching the minimum. 

We can illustrate this using a 
quadratic function of two vari- 
ables f{x) as shown in the fig- 
ure to the right. Beginning at 
some point Xq near the mini- 
mum and using a line-search to 
find the minimum of / in that 
direction to determine each step 
size (green line), an unfortunate 
starting position can result in 
a long zig-zag path and a large 
number of steps. 

To avoid this, we can use the conjugate-gradient 
method [36], which works by only stepping in directions 
that are conjugate to those already used. Writing a gen- 
eral quadratic function of many variables as 




Figure 21. Illustra- 
tion [50] of quadratic 
function (blue) mini- 
mization using the gra- 
dient descent (green) 
and conjugate gradient 
(red) methods. 



/(^) 



const. 



with X e C*, A e Mdxd, we can define A 
b — A^b such that: 



A^A and 



1 



'J Ax — x^b + const. 



A vector x is conjugate to another vector y with respect 
to / if and only if x^ Ay = 0. The gradient of the function 
is 

Vf{x) ^ Ax-b 

such that a stationary point x^ satisfies Ax^ — b. 
Given a basis consisting of n mutually conjugate vec- 
tors p\Apm — Snm, wc Can expand a;* in that basis 
Xt. = CnPn with coefficients dependent only on the 
corresponding basis vectors: 



lb 



plApn 



This means we can pick a starting vector po and proceed 
to the minimum in exactly n steps by finding successive 
Pn that are conjugate to all previous Po...n-i; where the 
step-size c„ is uniquely determined by the current direc- 
tion Pn- This can be further improved on by choosing 
specific Pn 



Pn+l 



PnPn 



where r„ = —\/f{xn) is the negative gradient and /?„ is 
a number defined below. We make steps 



Xn-\-l — Xn ~\~ OZ-nPn 



with 



-Vf{x 



n+l) 



CtnApn 



where po ^ ro — b — Axq- Requiring r}^r„ 
pl^Apjn = Snm then result in 



pliApn 



and Pn — 



SO that we only need to know p„ and r„ to calculate the 
next step. The red line in the above illustration demon- 
strates this procedure. We can modify this version of 
the conjugate-gradient method again so that it can itera- 
tively find the solution of approximately quadratic prob- 
lems. In this case, we must treat the function / and its 
gradient V/ as black boxes. We can do this by obtaining 
an approximate a, which is the only quantity requiring 
direct knowledge of A, by doing a line-search to find the 
minimum of f{xn + otpn)- The algorithm, known as the 
non-linear conjugate gradient method, is then: 

1. Calculate r„ — — V/(a?„). 

2. Compute Pn-i- 

3. Calculate the next conjugate vector p„ = r„ -|- 

Pn-lPn-l- 

4. Use a line-search to find a„ = argmiuQ, /(cc„ + 
apn). 

5. Set new position Xn+i = x.^ + cunPn- 

The initial values are Po = '''o = ~V/(a;o). If / is ex- 
actly quadratic and we ignore all numerical error, this 
algorithm will find the minimum in d iterations or less. 
With an approximately quadratic / and/or accounting 
for limited numerical precision, the vectors p„ will not 
be exactly conjugate to each other and errors will ac- 
cumulate. The algorithm must therefore be re-started at 
least every d iterations. Note that there are other choices 
for Pn that are equivalent in the quadratic case, but re- 
sult in different non-linear conjugate gradient algorithms. 
The above choice is the one originally used by Fletcher 
and Reeves [51]. 
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Appendix B: Matrix elements for phi-4-theory 



The matrix elements of the operators needed to imple- 
ment the (/)^-theory Hamiltonian using the site number- 
basis defined in section III are as follows: 



'2^2 



+ Ss,t+lt^/{t+l) 

+ Ss,t+i{{t + l) + {t + 2))^{t + l) 



+ 6s+iAis + 1) + (s + 2))V(s + l) 

-t- Ss+l^tS\/ (s + 1) 

+Ss+3,t V(s + l)(s + 2)(s + 3) 



-f(5s.t_iv/(s + l)l , 



+ SsA{Qt^ + Qt + 3) 



+ (5,+2,t(4s + 6)V(s + l)(s + 2) 
+'5s+3,t-i + + 2)(s + 3)(s + 4)1 , 



{s\n\t) 



V2 



'5.,t+iv/(* + l) 



-Ss+i,t\/ {s + 1) 



+ (5,,t(2s + l) 
+<5.+i,t-iV(-'^ + l)(s + 2) 



+ <5,,t(2s-)-l) 

-,5,+i,t_iV(s + l)(s + 2) . 



